#Librerias
library(IRdisplay)
library(ggplot2)
library(cowplot)
library(latex2exp)
library(data.table)
library(dplyr)
library(formattable)
library(tidyr)
library(afex)
| $n$ | $\sum_{i=1}^{n} x_{i1}$ | $\sum_{i=1}^{n} x_{i2}$ | $\sum_{i=1}^{n} y_{i}$ | $\sum_{i=1}^{n} y^2_{i}$ | $\sum_{i=1}^{n} x^2_{i1}$ | $\sum_{i=1}^{n} x^2_{i2}$ | $\sum_{i=1}^{n} x_{i1}x_{i2}$ | $\sum_{i=1}^{n} x_{i1}y_{i}$ | $\sum_{i=1}^{n} x_{i2}y_{i}$ |
|---|---|---|---|---|---|---|---|---|---|
| 10 | 223 | 553 | 1916 | 371595.6 | 5200 | 31729 | 12352 | 43550 | 104736.8 |
n <- 10
sum_x1 <- 223; sum_x2 <- 553
sum_y <- 1916; sum_y_sq <- 371595.6
sum_x1_sq <- 5200; sum_x2_sq <- 31729; sum_x1_x2 <- 12352
sum_x1_y <- 43550; sum_x2_y <- 104736.8
A <- as.matrix(cbind(
c(n, sum_x1, sum_x1_x2),
c(sum_x1, sum_x1_sq, sum_x1_x2),
c(sum_x2, sum_x1_x2, sum_x2_sq)))
B <- as.matrix(c(sum_y, sum_x1_y, sum_x2_y))
beta <- solve(A) %*% B
display_markdown(sprintf("$$\\begin{bmatrix} \\hat{\\beta_0} \\\\ \\hat{\\beta_1} \\\\ \\hat{\\beta_2}\\end{bmatrix} = \\begin{bmatrix} %f \\\\ %f \\\\ %f\\end{bmatrix}$$", beta[1], beta[2], beta[3]))
x <- as.matrix(c(1, 18, 43))
y <- sum(beta * x)
display_markdown(sprintf("$\\hat{y} = \\hat{\\beta_0} + \\hat{\\beta_1} x_1 + \\hat{\\beta_2} x_2 = %.4f$", y))
$\hat{y} = \hat{\beta_0} + \hat{\beta_1} x_1 + \hat{\beta_2} x_2 = 151.1852$
| $y$ | $x_1$ | $x_2$ | $x_3$ | $x_4$ |
|---|---|---|---|---|
| 240 | 25 | 24 | 91 | 100 |
| 236 | 31 | 21 | 90 | 95 |
| 290 | 45 | 24 | 88 | 110 |
| 274 | 60 | 25 | 87 | 88 |
| 301 | 65 | 25 | 91 | 94 |
| 316 | 72 | 26 | 94 | 99 |
| 300 | 80 | 25 | 87 | 97 |
| 296 | 84 | 25 | 86 | 96 |
| 267 | 75 | 24 | 88 | 110 |
| 276 | 60 | 25 | 91 | 105 |
| 288 | 50 | 25 | 90 | 100 |
| 261 | 38 | 23 | 89 | 98 |
df <- data.frame(read.table("./TP3_tables/data2.txt", header = TRUE))
names(df) <- c("y", "x1", "x2", "x3", "x4")
display_markdown("#### **Resumen de los datos**")
df_summary <- as.data.frame(apply(df, 2, summary))
df_summary$values <- rownames(df_summary)
df_summary <- df_summary[,c(6,1,2,3,4,5)]
colnames(df_summary)[1] <- " "
rownames(df_summary) <- c()
table <- formattable(df_summary, align='c', list(` ` = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL)
model <- lm(y ~ x1 + x2 + x3 + x4, data=df)
coef <- as.data.frame(x=model$coefficients)
colnames(coef) <- c('')
rownames(coef) <- c('β0', 'β1', 'β2', 'β3', 'β4')
coef <- as.data.frame(t(coef))
table <- formattable(coef, align = c("c", "c", "c", "c", "c"))
display_markdown("#### **Coeficientes**")
as.htmlwidget(table, width="50%", height=NULL)
x_test <- data.frame(x1=75, x2=24, x3=30, x4=98)
y_test <- predict(model, x_test, interval = "none", type = "response")
display_markdown(paste("$\\hat{y} = \\hat{\\beta}_1 x_1 + \\hat{\\beta}_2 x_2 + \\hat{\\beta}_3 x_3 + \\hat{\\beta}_4 x_4 = $", round(y_test, 4)))
$\hat{y} = \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \hat{\beta}_3 x_3 + \hat{\beta}_4 x_4 = $ 201.3144
f_stat<- as.data.frame(t(summary(model)$fstatistic))
f_stat$p_value <- c(pf(f_stat[,1], f_stat[,2], f_stat[,3], lower.tail=FALSE))
f_stat <- as.data.frame(cbind('F-statistic', f_stat))
colnames(f_stat) <- c(' ', 'value', 'df1', 'df2', 'p-value')
table <- formattable(f_stat, align='c', list(` ` = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL, )
y_pred <- predict(model) # valores predichos
y_true <- df['y'] # valores verdaderos
SSE <- sum((y_true - y_pred)^2)
n <- nrow(df) # n=12 observaciones
p <- ncol(coef) # p=5 coeficientes de regresión
MSE <- SSE / (n - p)
display_markdown(paste("$\\hat{\\sigma^2} = \\text{MSE} =$", round(MSE, 4)))
$\hat{\sigma^2} = \text{MSE} =$ 242.7156
# Otra forma
residuals <- anova(model)['Residuals',]
residuals$values <- rownames(residuals)
residuals <- residuals[,c(6,1,2,3)]
colnames(residuals)[1] <- " "
rownames(residuals) <- c()
table <- formattable(residuals, align='c', list(` ` = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL, )
model_coef <- format(round(model_summary$coefficients, 5), nsmall=5)
model_coef <- cbind(coef=c('β0', 'β1', 'β2', 'β3', 'β4'), model_coef)
colnames(model_coef)[1] <- "Coefficient"
rownames(model_coef) <- c()
table <- formattable(as.data.frame(model_coef, 5), align='c', list('Coefficient' = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL)
alpha <- 0.05
beta <- as.data.frame(model_summary$coefficients)['Estimate']
std_dev <- as.data.frame(model_summary$coefficients)['Std. Error']
t <- qt(alpha/2, 8, lower.tail=FALSE) # 7 grados de libertad
conf_int <- as.data.frame(c(round(beta - t * std_dev, 4), round(beta + t * std_dev, 4)), col.names=c('lwr', 'upr'), row.names=c('$\\beta_0$', '$\\beta_1$', '$\\beta_2$', '$\\beta_3$', '$\\beta_4$'))
conf_int <- conf_int[-1,]
display_markdown("Intervalos de confianza del 95% para los coeficientes de regresión:")
for (row in 1:nrow(conf_int))
{
display_markdown(glue::glue("{rownames(conf_int)[row]}: ({conf_int$lwr[row]}, {conf_int$upr[row]})"))
}
Intervalos de confianza del 95% para los coeficientes de regresión:
$\beta_1$: (-0.2453, 1.456)
$\beta_2$: (-3.2994, 21.1467)
$\beta_3$: (-4.0776, 6.9525)
$\beta_4$: (-1.6786, 1.7058)
names <- c('β1', 'β2', 'β3', 'β4')
conf_int$avg <- apply(conf_int, 1, mean)
# Gráfico de los intervalos de confianza
options(repr.plot.width=7, repr.plot.height=4)
ggplot(conf_int, aes(names, avg, colour=names)) +
ggtitle("Intervalos de confianza de 95% para los coeficientes de regresión") +
geom_errorbar(aes(ymin=lwr, ymax=upr), width = 0.2) +
labs(x=TeX('Coeficientes regresión ($\\beta$ )'), y='Intervalo de confianza', color='Coeficientes') +
geom_hline(yintercept=0, linetype="dashed", col="black") +
theme_light()
x_test <- data.frame(x1=75, x2=24, x3=90, x4=98)
y_conf <- as.data.frame(predict(model, x_test, interval="confidence", level=0.95, type="response"))
colnames(y_conf) <- c('ajuste', 'cota inferiror', 'cota superior')
y_conf <- select(y_conf, 2, 1, 3)
display_markdown('#### **Intervalo de confianza:**')
as.htmlwidget(formattable(as.data.frame(y_conf), align="c"), width="50%")
display_markdown(sprintf("Intervalo de confianza del 95%% para la media de la potencia mensual consumida: $\\left(%.4f, %.4f\\right)$", y_conf[1], y_conf[3]))
Intervalo de confianza del 95% para la media de la potencia mensual consumida: $\left(263.7879, 311.3357\right)$
y_pred <- as.data.frame(predict(model, x_test, interval="prediction", level=0.95, type="response"))
colnames(y_pred) <- c('ajuste', 'cota inferiror', 'cota superior')
y_pred <- select(y_pred, 2, 1, 3)
display_markdown('#### **Intervalo de predicción:**')
as.htmlwidget(formattable(as.data.frame(y_pred), align="c"), width="50%")
display_markdown(sprintf("Intervalo de predicción del 95%% para la media de la potencia mensual consumida: $\\left(%.4f, %.4f\\right)$", y_pred[1], y_pred[3]))
Intervalo de predicción del 95% para la media de la potencia mensual consumida: $\left(243.7175, 331.4062\right)$
r_squared <- summary(model)$r.squared
display_markdown(sprintf('$R^2 = %f$', r_squared))
$R^2 = 0.744750$
# LOWESS line (la misma linea que ajusta los valores usando plot(model)):
smoothed <- data.frame(with(df, lowess(x = model$fitted, y = model$residuals)))
# Gráficos
options(repr.plot.width=7, repr.plot.height=4)
res_vs_fit <- ggplot(model) +
geom_point(aes(x=model$fitted, y=model$residuals),color= '#ff9696', size=2) +
geom_path(data = smoothed, aes(x = x, y = y), col="#a399ff") +
geom_hline(linetype = 2, yintercept=0, alpha=0.2) +
ggtitle("Residuals vs Fitted") +
xlab("Fitted") +
ylab("Residuals") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))
qq_plot <- ggplot(model, aes(qqnorm(.stdresid)[[1]], .stdresid)) +
geom_point(na.rm = TRUE,color= '#ff9696') +
geom_abline(col="#a399ff") +
xlab("Theoretical Quantiles") +
ylab("Standardized Residuals") +
ggtitle("Normal Q-Q") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))
plot_grid(res_vs_fit, qq_plot, ncol = 2)
| $y$ | $x_1$ | $x_2$ |
|---|---|---|
| 193 | 1.6 | 851 |
| 230 | 15.5 | 816 |
| 172 | 22.0 | 1058 |
| 91 | 43.0 | 1201 |
| 113 | 33.0 | 1357 |
| 125 | 40.0 | 1115 |
df <- data.frame(read.table("./TP3_tables/data3.txt", header = TRUE))
names(df) <- c("y", "x1", "x2")
display_markdown("#### **Resumen de los datos**")
df_summary <- as.data.frame(apply(df, 2, summary))
df_summary$values <- rownames(df_summary)
df_summary <- df_summary[,c(4,1,2,3)]
colnames(df_summary)[1] <- " "
rownames(df_summary) <- c()
table <- formattable(df_summary, align='c', list(` ` = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL)
model_1 <- lm(y ~ x1 + x2, data=df)
display_markdown('#### **Primer modelo** $\\left(y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 \\right)$')
coef <- as.data.frame(x=model_1$coefficients)
colnames(coef) <- c('')
rownames(coef) <- c('β0', 'β1', 'β2')
coef <- as.data.frame(t(coef))
table <- formattable(coef, align = 'c')
as.htmlwidget(table, width="50%", height=NULL)
x_test <- data.frame(x1=25, x2=1000)
y_test_1 <- predict(model_1, x_test, interval = "none", type = "response")
display_markdown(paste("$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 \\times 25 + \\hat{\\beta}_2 \\times 1000 = $", round(y_test_1, 4)))
$\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 \times 25 + \hat{\beta}_2 \times 1000 = $ 165.2902
model_2 <- lm(y ~ x1 * x2, data=df)
display_markdown('#### **Segundo modelo** $\\left(y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 + \\beta_{1:2} x_1 x_2\\right)$')
coef <- as.data.frame(x=model_2$coefficients)
colnames(coef) <- c('')
rownames(coef) <- c('β0', 'β1', 'β2', 'β12')
coef <- as.data.frame(t(coef))
table <- formattable(coef, align = 'c')
as.htmlwidget(table, width="50%", height=NULL)
y_test_2 <- predict(model_2, x_test, interval = "none", type = "response")
display_markdown(paste("$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 \\times 25 + \\hat{\\beta}_2 \\times 1000 + \\hat{\\beta}_{1:2} \\times 25 \\times 1000 = $", round(y_test_2, 4)))
$\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 \times 25 + \hat{\beta}_2 \times 1000 + \hat{\beta}_{1:2} \times 25 \times 1000 = $ 184.4907
display_markdown('#### **Primer modelo** $\\left(y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 \\right)$')
f_stat<- as.data.frame(t(summary(model_1)$fstatistic))
f_stat$p_value <- c(pf(f_stat[,1], f_stat[,2], f_stat[,3], lower.tail=FALSE))
f_stat <- as.data.frame(cbind('F-statistic', f_stat))
colnames(f_stat) <- c(' ', 'value', 'df1', 'df2', 'p-value')
table <- formattable(f_stat, align='c', list(` ` = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL, )
display_markdown('#### **Segundo modelo** $\\left(y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 + \\beta_{1:2} x_1 x_2\\right)$')
f_stat<- as.data.frame(t(summary(model_2)$fstatistic))
f_stat$p_value <- c(pf(f_stat[,1], f_stat[,2], f_stat[,3], lower.tail=FALSE))
f_stat <- as.data.frame(cbind('F-statistic', f_stat))
colnames(f_stat) <- c(' ', 'value', 'df1', 'df2', 'p-value')
table <- formattable(f_stat, align='c', list(` ` = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL, )
display_markdown('#### **Primer modelo** $\\left(y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 \\right)$')
model_coef <- format(round(summary(model_1)$coefficients, 5), nsmall=5)
model_coef <- cbind(coef=c('β0', 'β1', 'β2'), model_coef)
colnames(model_coef)[1] <- "Coefficient"
rownames(model_coef) <- c()
table <- formattable(as.data.frame(model_coef, 5), align='c', list('Coefficient' = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL)
display_markdown('#### **Segundo modelo** $\\left(y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 + \\beta_{1:2} x_1 x_2\\right)$')
model_coef <- format(round(summary(model_2)$coefficients, 5), nsmall=5)
model_coef <- cbind(coef=c('β0', 'β1', 'β2', 'β12'), model_coef)
colnames(model_coef)[1] <- "Coefficient"
rownames(model_coef) <- c()
table <- formattable(as.data.frame(model_coef, 5), align='c', list('Coefficient' = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL)
full_model <- lm(y ~ x1 + x2, data=df)
sse_full <- sum((predict(full_model) - df$y)^2)
display_markdown(paste('$SSE_{MC} =$', round(sse_full, 4)))
$SSE_{MC} =$ 1950.4222
reduced_model <- lm(y ~ x1, data=df)
sse_reduced <- sum((predict(reduced_model) - df$y)^2)
display_markdown(paste('$SSE_{MR} =$', round(sse_reduced, 4)))
$SSE_{MR} =$ 3871.6309
# Cálculo de F0
k <- 2; r <- 1; n <- nrow(df); p <- 3
f0 <- ((sse_reduced - sse_full) / (k - r)) / (sse_full / (n - p))
display_markdown(paste('$F_0 =$', round(f0, 4)))
$F_0 =$ 2.9551
# Cálculo del p-valor
alpha <- 0.05
p_value <- 1 - pf(f0, df1=k-r, df2=n-p)
display_markdown(paste('$\\text{p-valor} =$', round(p_value, 4)))
$\text{p-valor} =$ 0.1841
display_markdown('Intervalos de confianza del $99\\%$ para $\\beta_1$ y $\\beta_2$:')
conf_int <- as.data.frame(confint(full_model, c('x1', 'x2'), level = 0.99))
conf_int <- cbind(c('β1', 'β2'), conf_int)
colnames(conf_int)[1] <- ' '
rownames(conf_int) <- c()
as.htmlwidget(formattable(conf_int, align='c', list(' ' = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left')))), width="50%")
Intervalos de confianza del $99\%$ para $\beta_1$ y $\beta_2$:
colnames(conf_int) <- c('lwr', 'upr')
conf_int$avg <- apply(conf_int, 1, mean)
coef <- c('β1', 'β2')
# Gráfico de los intervalos de confianza
options(repr.plot.width=7, repr.plot.height=4)
ggplot(conf_int, aes(coef, avg, colour=coef)) +
ggtitle("Intervalos de confianza de 99% para los coeficientes de regresión") +
geom_errorbar(aes(ymin=lwr, ymax=upr), width = 0.2) +
labs(x=TeX('Coeficientes regresión ($\\beta$ )'), y='Intervalo de confianza', color='Coeficientes') +
geom_hline(yintercept=0, linetype="dashed", col="black") +
scale_y_continuous(breaks=round(seq(min(conf_int) - 1, max(conf_int) + 1), 0)) +
theme_light()
inter_model <- lm(y ~ x1*x2, data=df)
display_markdown('Intervalos de confianza del $99\\%$ para $\\beta_1$, $\\beta_2$ y $\\beta_{1:2}$:')
conf_int <- as.data.frame(confint(inter_model, c('x1', 'x2', 'x1:x2'), level = 0.99))
conf_int <- cbind(c('β1', 'β2', 'β1:2'), conf_int)
colnames(conf_int)[1] <- ' '
rownames(conf_int) <- c()
as.htmlwidget(formattable(conf_int, align='c', list(' ' = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left')))), width="50%")
Intervalos de confianza del $99\%$ para $\beta_1$, $\beta_2$ y $\beta_{1:2}$:
colnames(conf_int) <- c('lwr', 'upr')
conf_int$avg <- apply(conf_int, 1, mean)
coef <- rownames(conf_int)
# Gráfico de los intervalos de confianza
options(repr.plot.width=7, repr.plot.height=4)
ggplot(conf_int, aes(coef, avg, colour=coef)) +
ggtitle("Intervalos de confianza de 99% para los coeficientes de regresión") +
geom_errorbar(aes(ymin=lwr, ymax=upr), width = 0.2) +
labs(x=TeX('Coeficientes regresión ($\\beta$ )'), y='Intervalo de confianza', color='Coeficientes') +
geom_hline(yintercept=0, linetype="dashed", col="black") +
scale_y_continuous(breaks=round(seq(min(conf_int) - 1, max(conf_int) + 1), -1)) +
theme_light()
r_squared <- summary(full_model)$r.squared
display_markdown(paste('$R^2 =$', round(r_squared, 4)))
$R^2 =$ 0.8618
r_squared <- summary(inter_model)$r.squared
display_markdown(paste('$R^2 =$', round(r_squared, 4)))
$R^2 =$ 0.9205
| $\beta_1$ | $\beta_2$ | $\beta_3$ | |
|---|---|---|---|
| $t_0$ | 4.82 | 8.21 | 0.98 |
n <- 25 # cantidad de observaciones
t_values <- c(4.82, 8.21, 0.98) # cocientes calculados
p <- length(t_values) + 1 # número de coeficientes = 4
df <- n - p # grados de libertad = número de coeficientes = 4
p_values <- pt(q=t_values, df=df, lower.tail=FALSE)
display_html(sprintf('
<table style="width: 40%%;">
<thead style="text-align: center">
<th style="text-align: center; font-weight: normal">Coeficiente</th>
<th style="text-align: center">$t_0$</th>
<th style="text-align: center">$\\text{p-valor}$</th>
</thead>
<tbody>
<tr>
<td style="text-align: center">$\\beta_1$</td>
<td style="text-align: center">%.4f</td>
<td style="text-align: center">%.e</td>
</tr>
<tr>
<td style="text-align: center">$\\beta_2$</td>
<td style="text-align: center">%.2f</td>
<td style="text-align: center">%.4f</td>
</tr>
<tr>
<td style="text-align: center">$\\beta_3$</td>
<td style="text-align: center">%.2f</td>
<td style="text-align: center">%.4f</td>
</tr>
</tbody>
</table>
', t_values[1], p_values[1], t_values[2], p_values[2], t_values[3], p_values[3]))
| Coeficiente | $t_0$ | $\text{p-valor}$ |
|---|---|---|
| $\beta_1$ | 4.8200 | 5e-05 |
| $\beta_2$ | 8.21 | 0.0000 |
| $\beta_3$ | 0.98 | 0.1691 |
| $y$ | 24.60 | 24.71 | 23.90 | 39.50 | 39.60 | 57.12 | 67.11 | 67.24 | 67.15 | 77.87 | 80.11 | 84.67 |
| $x$ | 4.00 | 4.00 | 4.00 | 5.00 | 5.00 | 6.00 | 6.50 | 6.50 | 6.75 | 7.00 | 7.10 | 7.30 |
x <- c(4.00, 4.00, 4.00, 5.00, 5.00, 6.00, 6.50, 6.50, 6.75, 7.00, 7.10, 7.30)
y <- c(24.60, 24.71, 23.90, 39.50, 39.60, 57.12, 67.11, 67.24, 67.15, 77.87, 80.11, 84.67)
data <- as.data.frame(cbind(x, y))
display_markdown("#### **Resumen de los datos**")
df_summary <- as.data.frame(apply(data, 2, summary))
df_summary$values <- rownames(df_summary)
df_summary <- df_summary[,c(3,1,2)]
colnames(df_summary)[1] <- " "
rownames(df_summary) <- c()
table <- formattable(df_summary, align='c', list(` ` = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL)
model <- lm(y ~ x + I(x^2), data=data)
coef <- as.data.frame(x=model$coefficients)
colnames(coef) <- c('')
rownames(coef) <- c('β0', 'β1', 'β11')
coef <- as.data.frame(t(coef))
table <- formattable(coef, align = 'c')
as.htmlwidget(table, width="50%", height=NULL)
# Gráfico del modelo
options(repr.plot.width=7, repr.plot.height=4)
ggplot(data, aes(x=x, y=y)) +
labs(
title="Modelo de regresión de segundo orden para Eficiencia del impulso",
subtitle=TeX("( $y = β_0 + β_1 x + β_{11} x^2$ )"),
x="Ángulo de divergencia [º]",
y="Eficiencia del impulso [%]") +
stat_smooth(aes(x=x, y=y, col="Curva de ajuste", fill="Intervalo de confianza (95%)"),
method="lm", formula=y ~ x + I(x^2), se=TRUE, size=1) +
geom_point(aes(shape='Datos del experimento'), col='#ff9696') +
theme_light() +
theme(
plot.title = element_text(size=10, hjust = 0.5),
plot.subtitle = element_text(size=10, hjust = 0.5),
legend.position = "bottom") +
scale_fill_manual(NULL, values = '#BAB5E3') +
scale_color_manual(NULL, values = '#a399ff') +
scale_shape_manual(NULL, values = 19) +
guides(
color=guide_legend(override.aes = list(fill=NA), order=1),
fill=guide_legend(override.aes = list(color=NA), order=2),
shape=guide_legend(order=3))
f_stat<- as.data.frame(t(summary(model)$fstatistic))
f_stat$p_value <- c(pf(f_stat[,1], f_stat[,2], f_stat[,3], lower.tail=FALSE))
f_stat <- as.data.frame(cbind('F-statistic', f_stat))
colnames(f_stat) <- c(' ', 'value', 'df1', 'df2', 'p-value')
table <- formattable(f_stat, align='c', list(` ` = formatter("span",style = ~ style(
'font-weight'='bold', 'text-align'='left'))))
as.htmlwidget(table, width="50%", height=NULL, )
r_squared <- summary(model)$r.squared
display_markdown(sprintf('$R^2 = %.4f$', r_squared))
$R^2 = 0.9957$
sse_full <- sum((predict(model) - data$y)^2)
display_markdown(paste('$SSE_{MC} =$', round(sse_full, 4)))
$SSE_{MC} =$ 24.7202
reduced_model <- lm(y ~ x, data=data)
sse_reduced <- sum((predict(reduced_model) - data$y)^2)
display_markdown(paste('$SSE_{MR} =$', round(sse_reduced, 4)))
$SSE_{MR} =$ 48.9825
# Cálculo de F0
k <- 2; r <- 1; n <- nrow(data); p <- length(coef)
f0 <- ((sse_reduced - sse_full) / (k - r)) / (sse_full / (n - p))
display_markdown(paste('$F_0 =$', round(f0, 4)))
$F_0 =$ 8.8333
# Cálculo del p-valor
alpha <- 0.05
p_value <- 1 - pf(f0, df1=k-r, df2=n-p)
display_markdown(paste('$\\text{p-valor} =$', round(p_value, 4)))
$\text{p-valor} =$ 0.0156
# LOWESS line (la misma linea que ajusta los valores usando plot(model)):
smoothed <- data.frame(with(data, lowess(x = model$fitted, y = model$residuals)))
# Gráficos
options(repr.plot.width=7, repr.plot.height=4)
res_vs_fit <- ggplot(model) +
geom_point(aes(x=model$fitted, y=model$residuals), color= '#ff9696', size=2) +
geom_path(data = smoothed, aes(x = x, y = y), col="#a399ff") +
geom_hline(linetype = 2, yintercept=0, alpha=0.2) +
ggtitle("Residuals vs Fitted") +
xlab("Fitted") +
ylab("Residuals") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))
qq_plot <- ggplot(model) +
stat_qq(aes(sample = .stdresid), color= '#ff9696') +
geom_abline(col="#a399ff") +
xlab("Theoretical Quantiles") +
ylab("Standardized Residuals") +
ggtitle("Normal Q-Q") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))
plot_grid(res_vs_fit, qq_plot, ncol = 2)
cubic_model <- lm(y ~ x + I(x^2) + I(x^3), data=data)
coef <- as.data.frame(x=cubic_model$coefficients)
colnames(coef) <- c('')
rownames(coef) <- c('β0', 'β1', 'β2', 'β3')
coef <- as.data.frame(t(coef))
table <- formattable(coef, align = 'c')
as.htmlwidget(table, width="50%", height=NULL)
# Gráfico del modelo
options(repr.plot.width=7, repr.plot.height=4)
ggplot(data, aes(x=x, y=y)) +
labs(
title="Modelo de regresión de tercer orden para Eficiencia del impulso",
subtitle=TeX("( $y = β_0 + β_1 x + β_2 x^2 + β_3 x^3$ )"),
x="Ángulo de divergencia [º]",
y="Eficiencia del impulso [%]") +
stat_smooth(aes(x=x, y=y, col="Curva de ajuste", fill="Intervalo de confianza (95%)"),
method="lm", formula=y ~ x + I(x^2) + I(x^3), se=TRUE, size=1) +
geom_point(aes(shape='Datos del experimento'), col='#ff9696') +
theme_light() +
theme(
plot.title = element_text(size=10, hjust = 0.5),
plot.subtitle = element_text(size=10, hjust = 0.5),
legend.position = "bottom") +
scale_fill_manual(NULL, values = '#BAB5E3') +
scale_color_manual(NULL, values = '#a399ff') +
scale_shape_manual(NULL, values = 19) +
guides(
color=guide_legend(override.aes = list(fill=NA), order=1),
fill=guide_legend(override.aes = list(color=NA), order=2),
shape=guide_legend(order=3))
sse_full <- sum((predict(cubic_model) - data$y)^2)
display_markdown(paste('$SSE_{MC} =$', round(sse_full, 4)))
$SSE_{MC} =$ 22.4152
reduced_model <- lm(y ~ x + I(x^2), data=data)
sse_reduced <- sum((predict(reduced_model) - data$y)^2)
display_markdown(paste('$SSE_{MR} =$', round(sse_reduced, 4)))
$SSE_{MR} =$ 24.7202
# Cálculo de F0
k <- 3; r <- 2; n <- nrow(data); p <- length(coef)
f0 <- ((sse_reduced - sse_full) / (k - r)) / (sse_full / (n - p))
display_markdown(paste('$F_0 =$', round(f0, 4)))
$F_0 =$ 0.8226
# Cálculo del p-valor
alpha <- 0.05
p_value <- 1 - pf(f0, df1=k-r, df2=n-p)
display_markdown(paste('$\\text{p-valor} =$', round(p_value, 4)))
$\text{p-valor} =$ 0.3909